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A COMPUTATIONAL SYSTEM FOR 
AERODYNAMIC DESIGN AND ANALYSIS OF 
SUPERSONIC AIRCRAFT 

PART 1 - GENERAL DESCRIPTION AND THEORETICAL DEVELOPMENT 

W. D. Middleton and J. L. Lundry 
Boeing Commercial Airplane Company 


1 . 0 SUMMARY 


An improved integrated system of computer programs has been 
developed for the design and analysis of supersonic 

configurations. The system uses linearized theory methods for the 
calculation of surface pressures and supersonic area rule concepts 
in combination with linearized theory for calculation of 
aerodynamic force coefficients. Interactive graphics are optional 
at the user's request. 

The integrated system consists of an executive "driver" and seven 
basic computer programs including a geometry input module, which 
are used to build up the force coefficients of a selected 
configuration. 

The description of the design and analysis system is covered in 
three separate documents with the following subtitles: 

Part 1 - General Description and Theoretical Development 
Part 2 - User's Manual 

Part 3 - Computer Program Description 

These documents supersede NASA Contractor Reports CR-2520, CR- 
2521, and CR-2522 which described an earlier version of the 
system. 

This document (part 1) presents a general description of the 
system and describes the theoretical methods used. 


1 




I 


2 . 0 INTRODUCTION 

Over a period of years, NASA-Langley has developed a basic 
computerized series of supersonic design and analysis methods for 
aerodynamic configuration studies (reference 1). The methods are 
characterized by their reliability in use and input simplicity. 

The Boeing Company has extended this basic series of methods and 
combined them into an integrated system of computer programs. The 
extensions to the methods provide several new features: 

• Addition of a near-field (thickness pressure) wave 
drag program, to complement the existing super.sonic 
area rule program. 

• Improved modeling of fuselage in lifting surface 
design and analysis programs. 

• Addition of configuration- dependent loadings in wing 
design program, so that the wing design is performed 
in the presence of fuselage and nacelle effects. 

• Addition of pressure limiting terms in the lifting 
pressure programs, to constrain the linear theory 
solution. 

• Optional CRT displays of selected program input and 
output data and provisions for limited user editing 
and intervention. 

A plot module is included in the system to produce configuration 
drawings, and a common geometry module is used to permit a single 
geometry input for all programs. 

The basis of the system is supersonic linearized theory, modified 
in two respects: 

• The "Whitham" correction to disturbance positioning 
is used in the propagation of body pressure fields. 

• The wing lifting pressure modules contain an 
optional limiting pressure feature to control the 
permissabie level of upper surface pressure 

coefficient. In addition, the wing design module contains 
a further constraint feature to limit the upper surface 
streamwise pressure gradient. 


Superposition is used to build up the theoretical force 
coefficients of a selected configuration. 
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The goals of the integrated system have been to develop an easily 
used supersonic design and analysis capability, with recognition 
of the need for constraints on linear theory methods to provide 
physical realism, and with inclusion of interactive display for 
increased design control over optimization cycles. 

The description of the design and analysis system is broken into 
three documents: 

Part 1 - General Description and Theoretical Development 
Part 2 - User's Manual 

Part 3 - Computer Program Description 

This report (part 1) presents a general description of the system 
and the theoretical methods used. 

These three documents supersede NASA contractor reports CR-2520, 
CR-2521, and CR-2522 which described an earlier version of the 
system. 
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3.0 SYMBOLS 


A (x) 


D 

f(y) 

h(Z) 

k 

L 

M 

N 

P 

q 

R 

R 


element fraction 
load strength factor 
cross sectional area at x 
trailing edge element fraction 
pitching moment reference length 
wing tip element fraction 
chord 

drag coefficient, D/qS 
skin friction coefficient 
lift coefficient, li'ft/qS 

pitching moment coefficient, pitching moment/qSc 
pitching moment coefficient at zero lift 
local pressure coefficient, (p - p )/q 

lifting pressure coefficient (lower surface pressure 
coefficient minus upper surface pressure coefficient) 

drag force 

Whitham function 

decay function for F ( y ) calculation 
Mach number parameter 

grid element variable in x ( length wise) direction 
Mach number 

grid element variable in spanwise direction 
pressure 

dynamic pressure, .7 M 2 
body radius 
influence function 

5 


I 



radial distance to point in flow field 

reference area or cross-sectional area (figure 4.3-3) 

temperature 

variable of integration (figure 4.3-4) 
free stream velocity 
x perturbation velocity 
body pressure decay function 
upwash 

lengthwise variable 

span station or distance variable 

vertical variable 

decay function parameter 

angle of attack 

local fuselage camber angle 

Wach number parameter = ^M 2 - 1 

x variable of integration 

surface shape (equation 13) or Lagrange multiplier 
(figure 4.4-4) 

span wise variable of integration 
density 

velocity potential 
radial angle variable 
region of influence 
dynamic viscosity coefficient 


free stream conditions 


incompressible condition (equation 5) 
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loading numbers (equation 31) 


i# j 

t tilted lengthwise variable (figure 4.3-5) 

Superscripts: 

» reference temperature condition (equation 2) or 

first derivative with respect to x 

„ second derivative with respect tc x 

# field point grid element 


7 




4.0 DISCUSSION 


The integrated supersonic design and analysis system is shown 
schematically in figure 4.0-1. It consists of an executive 
“driver" and seven basic computer programs. The individual 
programs, or modules, provide data for configuration design or 
analysis as follows: 

• Skin friction is computed using turbulent flat plate 
theory. 

• Wave drag is calculated from either far-field 
(supersonic area rule) or near- field (surface 
pressure integration) methods. The far-field 
method is used for wave drag coefficient calculations 
and for fuselage optimization according to area 

rule concepts. The near-field method is used 
primarily as an analysis tool, where detailed 
pressure distributions are of interest. 

• Lifting pressures, drag-due-to-li ft , pitching moment 
and trim drag are computed from the lift analysis 
program, which breaks arbitrary wing/fuselage/ 
canard/nacelles/horizontal tail configurations into a 
mosaic of "Macn-box" rectilinear elements which are 
employed in linear theory solutions. A complementary 
wing design and optimization program solves for the 
wing shape required to support an optimized pressure 
distribution at a specified flight condition. 

• The plot module draws configuration pictures according 
to size and view inputs. 

• The geometry module handles configuration geometry 
for the system. 

• Interactive graphics, for data display or editing, 
are optional in the system at user's request. 

Operating in the analysis mode, the force coefficients of a 
selected configuration are built up through superposition as shown 
in figure 4.0-2. 

The entire design and analysis system is a single overlaid 
program. The executive level of the system controls module 
execution by means of special identification cards in the input 
data. Transfer of data in the system between modules is handled 
by disk storage and common blocks. 

All configuration geometry is read (or updated) in the geometry 
module. Geometry inputs are patterned after those of the NASA-LRC 
configuration plotting program. All "paneling" of the 


9 




Drag-due-to-lift 

Zero Lift Wave Drag (Wave and vortex) 


FIGURE 4.0- 1. -INTEGRA TED SUPERSONIC DESIGN AND ANALYSIS SYSTEM 






SUPERPOSITION METHOD OF DRAG ANALYSIS 


Drag 



ZERO LIFT WAVE DRAG 


FIGURE 4.0-2.~ DRAG BUILDUP 
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configuration for theoretical analyses is accomplished within the 
programs, and the user prepares only "drawing- type" geometry. 

Additional data on system structure and input and output formats 
are presented in parts 2 and 3. The theoretical methods used in 
the modules are described in the remainder of this document. 


4.1 Skin Friction Drag 


The skin friction drag of an arbitrary configuration is calculated 
from the T method of reference 2. Smooth flat plate, adiabatic- 
wall, and turbulent boundary layer conditions are assumed, with 
transition assumed to occur at the leading edge of each 
conf iguration component. 

The theory and experimental verification of the T f method are 
given in detail in reference 2. The essential elements of the 
method are presented in Appendix C of the reference, and are 
summarized here for completeness. 

The T 1 method is based on the calculation of a compressible skin 
friction coefficient, C F , from a reference skin friction 
coefficient, c’p , for a selected 'Mach number, M*,, Reynolds number, 
Roor and adiabatic wall temperature, T w . (The subscript «> is used 
throughout to denote free stream conditions) . 

The wall temperature ratio, Tw/T^, , is calculated from one- 
dimensional relationships assuming a wall recovery factor of .89. 


«r- = 1 + 0.178 


(1) 


The prime superscript (T* , R T , etc.,) refers to conditions at 
which incompressible flow relations must be evaluated in order to 
represent compressible flow. Sommer and Short in reference 2 
obtained the key relationship 

■f" = 1 + 0.035 Mj + 0.1(5 Q*-- l) (2) 

^oo x oo 


The Reynolds number relationship is 


R' 1 



( 3 ) 


where the viscosity ratio is given by the Sutherland equation 
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with T in degrees Rankine. 


The incompressible skin friction coefficient is given by the 
Karman-Schoenherr equation 


which gives 
analogy: 


0.2b2 

v'cIT 


= log 
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the corresponding relationship 

0.2U2 


vc; 


= log 


10 


(V *’) 


( 5 ) 


based on the T’ 


( 6 ) 


This equation is solved iteratively for C p ' . The desired 
compressible skin friction coefficient is obtained from C ' by 


C 


F 




(T) 


Most of the skin friction program is involved with computing 
wetted areas and reference lengths for the configuration. 
Components which may exhibit large variations in reference length 
(such as the wing, tail, canard etc.,) are broken into strips in 
order to calculate C p accurately. 


4.2 Far-Field Wave Drag 


The far-field wave drag program computes the zero lift wave drag 
of an arbitrary configuration by means of the supersonic area 
rule. A description of the program has been given in reference 3, 
and is summarized here. The method assumes not only an integrated 
total lift of zero, but local lift everywhere zero, and thus 
includes none of the wave drag contribution associated with the 
generation of lift. 

The supersonic area rule is a generalization of the transonic area 
rule, which states that the transonic wave drag of an aircraft is 
essentially the same as the wave drag of an equivalent body of 
revolution having the same cross-sectional area distribution. In 
the supersonic area rule procedure, several equivalent bodies of 
revolution are produced by passing a series of parallel cutting 
planes through the configuration, as shown below. The cutting 


planes, are inclined with respect to the aircraft axis at the Mach 
angle m , and a single equivalent body is produced for the series 
of cutting planes at a constant azimuthal angle, e . The area of 
the equivalent body at each station is the projection of the area 
intercepted by the cutting plane onto a plane normal to the 
aircraft axis. 



The wave drag of each equivalent body is calculated from the von 
Karman slender-body formula (reference 4) , which gives the drag as 
a function of the free stream conditions and equivalent-body area 
distribution. 


D(e) 



A" (x.^ A” (x 2 ) log 




( 8 ) 


where and X 2 are lengthwise variables of 
integration and A M is the second derivative of the 
body area distribution. 


The wave drag of the aircraft at a given Mach number is calculated 
from the integrated average of the equivalent -body wave drags. 




2tt 

D(e) ae 


(9) 


A useful feature of the supersonic area rule occurs in the 
optimization of a fuselage area distribution to minimize the wave 
drag of a wing-f us elage combination. According to Sheppard 
(reference 5), the wave drag of a wing- fuselage can be written as; 



wing 



wing eq. 
+ body 


-(f) 

body wing eq. 

body 


( 10 ) 
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where the wing equivalent body is a body of revolution obtained 
by averaging from 0 = 0 to 2 n all the projected areas intercepted 
by cutting planes passing through the wing for each X station. 

Considering a wing-f uselage, then, and assuming that the wing 
geometry is fixed, both the first and third terms of equation 10 
are fixed. This leaves optimization of the second term only, so 
that to design the minimum drag configuration, the fuselage must 
be contoured to produce a wing- equivalent-body-plus-fuselage 
having the shape of a minimum drag body of revolution. 

The wave drag program uses the Eminton-Lord fairing through a 
discrete set of points, which defines the shape of a minimum drag 
body through the points. Using this fairing, and identifying a 
few '’control-points*' on the fuselage area definition, an optimized 
fuselage area definition may be obtained as the difference between 
the wing equivalent body and the combined wing-equivalent -body- 
plus- fuselage. 

The wave drag program has mechanized this solution, and also 
includes the effects of nacelles, tail, etc. , which are included 
similarly. The shape of the minimum drag fuselage area 
distribution, and the drag of the complete configuration including 
the minimum drag fuselage, may be calculated for each 
configuration input. 

Because of the configuration generality that can be handled, the 
far-field program is the primary source of zero-lift wave drag 
coefficients in the design and analysis system. A complementary 
near-field program, more restricted in scope but which computes 
detailed thickness pressures, is descrioed in the following 
section. 


4.3 Near-Field Wave Drag 


Isolated Component Pressures 


The near-field wave drag program computes zero-lift thickness 
pressure distributions for wing- body- nacelle-empennage 

configurations in supersonic flow. The pressure distributions are 
integrated over the cross-sectional areas of the configuration to 
obtain the resultant drag force. Three basic calculations are 
performed to obtain the required pressure fields; 

• Thickness pressure distribution for a wing or tail 
of arbitrary section and planform. 

• Thickness pressure distribution on surface of 
fuselage or nacelles. 
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• "Whitham” near-field calculations to define pres- 
sure distributions propagating from fuselage or 
nacelles. 


Superposition is used to calculate the interference drag terms 
associated with the pressure field of a given component acting on 
the surfaces of the other components. 


The surface pressure distributions calculated include some (from 
the nacelles) which contribute lift due to a non -symmetrical 
distribution of volume, so that the drag calculations which 
include these loadings contain wave drag due to lift 
contributions. However, by selection of program options and 
subtraction of certain drag contributions, a zero lift wave drag 
comparable to the far-field wave drag value may be obtained, as 
described on page 3 6 . 


Wing or em pennage thickness pressures. - The surface pressure 

coefficient on the upper (or lower) surface of a flat-mean-line 
wing of symmetrical surface shape is obtained by first calculating 
the corresponding velocity potential, differentiating with respect 
to x (to get u) , and then computing the pressure coefficient from 
the linear theory approximation, Cp = -2 u/Uqo. 


The velocity potential computation, from reference 13, is: 


where 


<f> (x, y) 
X 



<t> U>y) = - 


1_ ff X dn dE 

* j) V U-t) 2 - 8 2 (y-n) 2 

T 


= velocity potential at a defined wing field 
point (x, y) 

= surface slope (dz/dx) of wing section at 
a wing integration point 
- V»2 _ 1 

= x variable of integration 
= y variable of integration 

= subscript denoting interval of integration 
(surface of the wing planform within the 
Mach forecone from x,y) 


( 11 ) 


The wing thickness pressure coefficient is 


where 


C p (x,y) 




u _ p d< P (*»y) 

U 3x 


? = local pressure at x,y 

= free stream static pressure 


( 12 ) 
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u = x perturbation velocity 
u oo,q» = free stream velocity and dynamic pressure 

Integration of the velocity potential equation is performed by 
representing the wing as a grid of rectangular elements, and 
substituting a numerical summation for the integration. The wing 
element "Mach-box" system exactly follows the arrangement 
described in Section 4.4 (for the lifting pressure case) and is 
summarized in figure 4.3-1. The grid elements, identified by L 
and N, are defined such that L is numerically equal to x and N is 
numerically equal to 0y, where x and 0 y take on only integer 
values. Partial grid elements along the wing leading and trailing 
edges are used to permit a closer approximation to the actual wing 
planform. The grid system of figure 4.3-1 is rather coarse for 
illustrative purposes; in actual usage, many more grid elements 
are employed. 

With respect to a specified field point x,y, the upstream region 
of influence t is approximated by the shaded grid elements in 
figure 4.3-1. Each of the shaded elements has associated with it 
an influence factor, R, which relates the effect of the element 
and its surface slope to the total velocity potential at the field 
point. 


The factor R is obtained by an approximate solution to Equation 11 
as follows: 


Let y-r? = u so that du = -drj 

4> (x,y) = 


vm 


du d( 


(x -£) 2 - 3 2 u 2 


Substituting an average value of (x -£) = (x -£) and integrating 

with respect to u. 


Restoring y - 1 = u 


r -i 

, x 1 X . -1 Bu ! 

* (x,y) = 7 [ T sin 7^Tj u 

j 

* (x>y) x sin_1 i^fr] 


For a given grid element in the Mach forecone, the integration 
interval to n 2 gives 

0(y- n, ) = N* - N + 0.5 
0 (y “ ^2 ) = N* - N -0.5 
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L + 0.5 


Also, X-* = L* - 

L* and N* define the field point grid element and L and N define 
an element in the forecone- The numerical summation then becomes 


* (x * y) = TbE x l,n 



N* - N - 0.5 
L* - L + 0.5 


. -1 
sin 


N» - N + 0.5 ~| 
L* - L + 0.5J 


(1 6) 


with R being the quantity between the brackets for the element 

(L,N) . 


The character of the R function is such that the sum of R for all 
elements across the Mach forecone at a constant L value is equal 
to - ir . In considering the case of an infinite rectangular wing 
of constant thickness shape (wedge section) , it is seen that the 
- it summation must be the case- Since the <P function for 
successive L* rows is increased by 


_X 

TT0 


[IR] 


L=const . 


so that 


Ax 


1 

*3 


(-*)=- 


A 

3 


and 



2 


A 

6 


which is the established linear theory result. Values of R within 
the Mach forecone from a typical field element are presented in 
f igure 4.3-2. 


The computer program for the wing thickness pressure calculations 
sets up the wing grid system, defines the surface slope A for 
each grid element, performs the summation for 0 at selected field 
elements, and differentiates 0 numerically to obtain Cp . Wave 
drag is obtained by summing the pressure coefficients times the 
corresponding surface slopes of the wing. 


In the wing surface slope definition, slopes of the elements 
occupying space covered by the fuselage (if there is one) are 
set equal to zero to eliminate their drag contribution. Partial 
elements are used at the wing-body intersection to improve the 
geometry definition. 


As was found to be the case in lifting surface programs involving 
the grid system used in the wing calculations, some oscillation in 
computed pressure coefficient values occurs. This has been 
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suppressed by a smoothing of the calculated pressure coefficient 
values, using a 5 term equation of the form; 


.333 A T# _ C + .66? A_ # C + A 

L *- 2 P L »-2 L 1 P L»-1 "L* "L..1 

p .333 + *667 Ay,_^ + Ap, + 1.0 


(17) 


where A = factor defining size of associated element (1 if whole 
element, fraction if not) 


Fuselage and nacelle thickness pressures. - The pressure 

distribution on the surface of fuselage or nacelles is obtained 
from a method developed at the Boeing Company by R. M. Kulfan, 
based on the Lighthill theory of reference 6. The method is 
applicable to bodies having either smooth area distributions or 
bodies with slope discontinuities. As used in the near- field 
program, smooth area distributions are assumed except {if 
required) at the nose or aft end of the body. Open-nose bodies, 
such as nacelles, are permissable. 


The solution technique is summarized in figure 4.3-3. An axial 
perturbation velocity is calculated which is a function of the 
body cross- sectional area growth (and radius distribution) and a 
decay function which relates area growth to its effect on a given 
field point. 


The axial perturbation velocity, u, is given by 
u(x) = - 2 -i | X Mil a S' 


where 


3 R 


(O 


(O 


( 18 ) 


x = body field station 
W (Z) = decay function 

£ = x variable of integration 
Z - position function = x ~ £ 

0R , 

= body radius at $ 

S'? = first derivative of body cross-sectional 

** area (S) at £ 


‘(€) 


The pressure coefficient is calculated from one-dimensional flow 
relationships (rather than the simpler linear theory 
approximation), as shown in figure 4.3-3. The decay function, 
W (Z) , is tabulated on page 23 . 


Fuselage and nacelle flow-fie ld pressures. - The pressure field 
propagated by the fuselage or nacelles is calculated using the 
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K) 

1-0 


p Integration 
i interval 



PERTURBATION 

VELOCITY 



S = cross-sectional area 


where: 

dS may include dS smoo th aR d dS discontinuous 


UlX) = Au (xUooth 


+ £ Au 


(x> 


disc 


u(x) 


±J 


Au ; 


smooth 2 it 

AR'(xi) 


W(Z) 

0 ^ R <£) 


dS' 

dx 


d£ 


(x) 


disc 


AR'(xi) 

' . 0 . 


W(Z: 


Where integration 
excludes ds'^j^ 



j 


r 

f \ " 

3.5 I 

PRESSURE 

COEFFICIENT 

1 

P W 0.7 j 

! 

1.0 + 0.2 | 

(l- 0+R '(x)1(l-0 + u, x ,) 2 -1.0j 

-1.0> 


FIGURE 4.3-3.-BODY SURFACE PRESSURES (LIGHTHILL THEORY) 






Whitham solution of reference 7. This solution, which is 
converted to a method amenable to the digital computer in 
reference 8, calculates shock wave positions and strengths through 
a modification of linear theory re suits - 


The Whitham solution begins with the calculation of a 
F (y) , which is dependent on body geometry: 


where y 

h (Z) 


S (t) 


end 




1 _ . 

2 

h (Z) d S 


(t) 


body field station 

decay function (simi lar to body pressure 
calculation) 

y variable of integration 
t position function 

first derivative of body cross-sectional 
area (S) at t 

y at end of integration (see figure 4.3-4) 


f unc ti on. 


(19) 


The solution for F (y) is very much like that of the body thickness 
pressure solution, with the exception that the integration is 
carried out to the point (y enj j ) at which the aft-running Mach line 
from the body centerline at y emerges from the body contour. The 
F(y) calculation is summarized in figure 4.3-4, and the decay 
function h(Z) is tabulated on page 2 T • 


Development of the near-field pressure signature at a radial 
distance (r) from the body centerline is summarized in figure 4.3- 
5. The F (y) function is tilted, as shown, by displacing the F(y) 
function according to 


where 


Ay = k/F F(y) 


k = 


_ 1.6973 M 
B 1 ’ 5 


( 20 ) 


This process results in a double-valued pressure signature where 
shock waves are present, a situation which is corrected by an 
area-balancing technique. The area-balancing (shaded lobes are 
equal in area) defines shock wave locations and strengths. The 
resultant tilted/area-balanced signature is then converted to 
pressure coefficient by 


C 

P 


P ~ Poo = 2 F(y) 
^co V 2$r 


( 21 ) 
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DECAY FUNCTION 



F (y) = 


27T 




h (Z) dS' (t) 


S = Cross-sectional area 


PRESSURE 

COEFFICIENT 


c = p - p °o _2 F]y) 
P q °° Y^r 


where: 

r = radial distance from body 
F(y) function displaced, shock waves 
located as shown in fig. 4.3-5 

FIGURE 4.3-4.- WHITHAM F(y) FUNCTION (NEAR-FIELD PRESSURE BASIS ) 


to 
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The position of the pressure signature at radius (r) is then given 
by 

x = y + Br (22) 

where x is the longitudinal distance from the body nose. The 
"tilting' 1 thus produces a correction in pressure signature 
positioning, relative to simple translation along a Mach line, 
which is the result of the remarkable theory of reference 7. 


Interference Calculations 


The previous sections summarized the basic pressure calculations 
performed by the near-field wave drag program. Integration of the 
wing and fuselage or nacelle pressure distributions over their 
respective surfaces gives the "isolated" wave drags of the 
components. Interference drags are obtained by superposition; 
i.e., the pressure field of each component is imposed upon the 
surfaces of the other components to calculate the resultant 
interference forces. Interference terms between empennage 
surfaces and the rest of the configuration have been considered 
small, however, and are neglected. 


The computer solution allows for up to 3 pairs of nacelles located 
external to the wing-f uselage (or 2 pairs plus a single nacelle at 
Y=0) . The nacelles may be either above or below the wing (or 
both) . 


Wing- fuselage interference. - Fuselage -on-wing interference is 

obtained by calculating the near- field pressure signature from the 
fuselage at selected spanwise stations, and imposing them upon the 
corresponding wing sections as a bouyancy field. The spanwise 
stations are the same as those used in the isolated wing thickness 
pressure calculations. 


Wing-on- fuselage interference is obtained by computing wing 

thickness pressures in the area occupied by the fuselage, after 
the wing surface slopes in the fuselage area have been set equal 
to zero. The thickness pressures thus calculated are "carry-over" 
pressures, and are imposed upon the fuselage surface slopes by 
transposing them aft along Mach lines to the fuselage. The 
fuselage area covered by the wing root is deleted from the 
interference term. 
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Typical results from wing-body calculations are presented in 
figure 4.3-6. The wing data, both isolated and fuselage-on-wing 
interference, are converted to section drag form. 


Nacelle interference terms are calculated like the wing-body 
interference terms, with the provision that two alternate 
solutions may be specified for the nacelle-on-wing interference 
term. Available experimental data do not make it clear whether a 
"wrap” or “glance” solution, as shown in figure 4.3-7, is more 
correct. Since the nacelle-on -wing interference term is 
substantial, both solutions are available in the program 
(controlled by an input code) . 

The nacelle interference terms are calculated as follows: 

Na ce lle-on-win g te rm. - The nacelle pressure field acting on the 
wing is obtained by calculating nacelle pressure signatures at the 
same spanwise stations used for the wing thickness pressures (plus 
extra stations immediately adjacent to the nacelle centerlines) , 
then defining a composite signature by summing together the 
contributions from all nacelles. The nacelle pressure 
coefficients are doubled to account for reflection from the wing 
surf ace. 

In the case of nacelles both above and below the wing, separate 
nacelle pressure fields for nacelles below, and then above, the 
wing are calculated and then summed. 

A unique feature of the near- field approach lies in the solution 
choices available, since the pressure signatures generated by a 
nacelle to act on the wing surface may be terminated on 
encountering another nacelle or may be allowed to pass 
undiminished around (or through) it. In the first case, the 
"glance” solution, the nacelle pressure field is assumed to 
reflect from other nacelles in such a way that, the reflected field 
exerts no further influence upon the wing surface. In the second 
case, the "wrap" solution, the nacelle pressure field is assumed 
to propagate around the other nacelles as if no obstruction were 
offered, i.e., as if they were transparent to the pressure field. 

A comparison of the theoretical nacelle pressure field with 
experimental data is presented in figures 4.3-8 and 4.3-9 at Mach 
numbers of 2.7 and 1.1. At M =2.7, there is no difference 
between wrap and glance solutions. At M = 1.1, the solution 
choice appears to favor the glance approach except on the outboard 
wing. However, the experimental pressure points are too sparse to 
permit firm conclusions. 

Nacelle- on- fuselage term. - The effect of the nacelles on the 

fuselage area distribution is built up by integrating each 
nacelle’s pressure signature upon the fuselage area growth, with 
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PRESSURES “GLANCE" AWAY FROM WING AT ADJACENT 
NACELLES 



The nacelle pressure field and accompanying shock wave? "glance” 
away from the wing when encountering adjacent nacelles. In 
application, the nacelle generated pressure field is terminated on 
encountering another nacelle. 



PRESSURES "WRAP" AROUND ADJACENT NACELLE 

The nacelle pressure fields and accompanying shock waves 
"wrap" around adjacent nacelles. In application, the nacelle 
generated pressure field is allowed to pass through another 
nacelle as if it were transparent. 


FIGURE 4.3-7.-NACELLE GENERATED WING PRESSURE FIELD CONCEPTS 
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the provision that nacelles below the wing affect only the 
fuselage area below the wing in the region of the wing (and vice- 
versa for above-wing nacelles). In the fuselage region below the 
wing, or above the wing, the nacelle pressures are doubled to 
account for reflection (as rs done in the nacelle pressure field 
case) . 


^ysslage-on^nacelle term. - The fuselage pressure signature is 

imposed cn each nacelle surface as a bouyancy force. 


" Wing thickness pressures at span stations 
where nacelles are located are transferred aft along Mach lines 
from the wing surface to the nacelle centerline to obtain the 
buoyancy field acting on each nacelle. 


Nacelle-on-nacel le term. - The interference term of other nacelles 
acting on a selected nacelle is calculated by building up the 
composite buoyancy field, and then imposing it upon the nacelle 
surface. In the case of nacelles on opposite sides of the wing, 
the pressure signatures are cut off where intercepted by the wing. 

The nacelle-generated pressure field is assumed to pass 
undiminished around or through other nacelles that may be present 
when a particular nacelle pair interference term is being 
calculated much as in the "wrap” solution employed for nacelle on 
wing terms. The "wrap" or "glance" option is not provided in this 
case because of the generally negligible difference in results. 

Nacelle. im age ef fe cts. - If the nacelle is located next to the 
wing (or body) , an "image" nacelle is used to create the reflected 
buoyancy field. The reflected field is cut off forward or aft of 
the wing if the nacelle is not entirely under the wing. 

The principal image effect is caused by the nacelle pressure field 
reflecting off the wing back onto itself. However, the reflected 
nacelle pressure field acting on other nacelles is also computed. 
This solution also utilizes an image nacelle representation, and 
the reflected signature is restricted to whatever part of the 
"real" signature encountered a reflecting surface. 


Composite Thickness Pressure Signatures 
and Configuration Drag 

Typical theoretical thickness pressure distributions acting on the 
conf iguration components are presented in figure 4.3-10. The net 
effect of these pressure signatures on drag is obtained by summing 
the various inputs and integrating them over the surface of the 
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fIGURE 4.3-10.- NEAR-FIELD PRESSURE SIGNATURES, 1.1 "WRAP" SOLUTION 
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configuration. The corresponding drag terms are summarized in the 
data of figure 4.3-11. 

The near-field wave drag program has three calculation features 
that have no direct counterpart in the far-field wave drag 
program. These are; 

1) wrap or glance solution for nacelle pressure field acting on 

wing. The wrap solution is nominally the same as the far- 

field solution. 

2) Nacelle image effects. Because of the transparency 

assumption of the far-field program pressure propagation, 
nacelle pressures do not reflect off adjacent surfaces. 

(Through addition and subtraction of separate calculations 
involving "image” nacelles, however, comparable drag data can 
be generated using the far-field program. ) 

3) Above- or below- wing fuselage area separation in the 
nacelle-on-fuselage term. Directly comparable results are 
obtained only for mid-wing configurations. 

Using the "wrap” solution, a mid-wing arrangement, and subtracting 
the nacelle image drag terms, a direct comparison between near- 
field and far-field program drag calculations can be made. 


4.4 Drag-Due-to-Lift (Design and Analysis) 

The wing design and lift analysis modules are separate lifting 
surface methods which solve the direct or inverse problems of; 

• Design - to define the wing camber surface shape 
required to produce a selected lifting pressure 
di stribution. The wing design program includes methods 
for defining an optimum (least drag) pressure 
distribution. 

• Lift analysis - to define the lifting pressure 
distribution acting on a given wing camber surface, and 
calculate the associated force coefficients. 


Originally, these programs were developed to utilize the "Mach- 
box" method of wing representation and evaluation of linear theory 
integral equations (references 9 through 11). In the Mach - box 
method, the wing is replaced by a grid of small rectilinear 
elements. Since many elements (in the thousands) can be used, a 
detailed description of complicated surface shapes, with 
associated computational accuracy, is possible. 
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The wing design and lift analysis programs have been expanded 
several times to add features (e.g., reference ll) and were 
recently reviewed in reference 12. The discussion of the 
aerodynamic theory of the Mach-box method in reference 12 is quite 
detailed, and is summarized here for completeness. The numerical 
method for the '’direct” case of wing design and optimization is 
given first, followed by the inverse case of lift analysis on page 
5b . 


Design Case 


Camber surface for a given loading. - A. typical wing planform 
described by a rectangular cartesian coordinate system is 
illustrated in figure 4.4-1. The y co ordinat e has been multiplied 
by the Mach number parameter B = Y m 2 - l for convenience in 
mathematical manipulations. 


For a wing of zero thickness lying essentially in the z~0 plane, 
linearized theory for supersonic flow defines the wing surface 
shape necessary to support a specified lift distribution by the 
inteqral equation 

(23) 


3z 
c 

3x 


(x,y) 
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iC P (x,y) + "fc 



(x - o AC u,n) dn 
E 

(y - h ) 2 l/(x - C ) 2 - $ 2 (y - n ) 2 


which is slightly modified form of equation (77a) of reference 13. 
With respect to a specified field point x, y, the upstream region 
of influence, t , enclosed by the Mach forecone is shown by the 
shaded area. 


Equation 23 may be rewritten into an influence function form by 
introducing the factor P, such that 


3z 

3^(x,y) 


-6 

~T 


AC p (X’y> 



R(x-s>y-n) ACp(5,n) d&n 


(2U) 


where R is defined as 

R(x-£»y-n) » * * ■ — ^ — — - 

e 2 (y - n ) 2 V (x - c ) 2 - e 2 (y - n ) 2 


(25) 


R is a function relating the local loading at point £ , to its^ 

influence in determining the downwash (or upwash) at the field 
point. 
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For the numerical evaluation of this integral, the wing is 
replaced by the grid system of figure 4,4- 2 • The grid elements, 
identified by L and N, are defined such that L is equal to X and 
N is equal to 0y , where X and 0y take on only integer values. 
Partial grid elements along the wing leading and trailing edges 
are used to permit a closer approximation to the actual wing 
planform- The grid system of figure 4.4-2 is rather coarse for 
illustrative purposes. In actual usage, many grid elements are 
employed. 


In the grid system, the field point x, y is located at the aft 
midpoint of a field point element L*, N* , and the region of 
integration, t , is approximated by the shaded grid elements. 


The integration required by equation 24 is performed for each 
element within the Mach forecone, considering Ac p constant over 
the element and using an average value of X-£ - L*-L+5. The 
resulting expression, derived in reference 12, is the factor R 
which is the value of the influence function for an element L, N: 


5, T „„ V(L« - L + 0.5) 2 - (K» - H - 0.5) 2 

R(L -L,N -N) - (l* - L + 0 .5)(N* - N - 0.5) 


_ V 


(L* - L + 0.5 ) 2 - (N* - N + 0.5 ) 2 
(L* - L +0.5) (N* - N + 0.5) 


(26) 


The integral equation 23 may then be replaced by the summation 
given below, where the summation includes all elements within the 
forecone and on the wing planform. The factors A, B, and C are 
element fractions for the wing leading edge, trailing edge, and 
wing tip, respectively, to allow for partial elements at those 
locations. 

= 4- AC 

ax * p ( 27 ) 

g R(L*-L,N*-N> A(L,N) B(L,N) C(L,N) AC (L,N) 

+ Ti p 

T 

The character of the R function is such that, f or_ a given L*-L set 
of elements within the forecone, the sum of the R values is zero, 
the single negative value at N*-N =0 balancing all the others (see 
figure 4.4-3). _At L* =L, where there is only one element in the 
summation, the R value is zero. 


The physical significance of this R variation is that (for 
positive lift) , all elements directly ahead of the field point 
element contribute only downwash and all other elements contribute 
upwash. An element at the leading edge near the wing tip of a 
subsonic leading edge wing, therefore, sees a concentrated upwash 
field. It is this upwash field that makes the subsonic leading 
edge twisted and cambered wing attractive from the standpoint of 
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drag-due- to- lift , since a local element: may be inclined forward to 
produce both lift and thrust. Conversely, an element located at 
the trailing edge and centerline of a diamond planform produces 
lift inefficiently, since the element sees an influencing region 
that contributes downwash. 


Ih computing the local surface slope required to obtain a 
specified lifting pressure coefficient, the ACp value is known 
everywhere on the wing from loading formulas. In the wing design 
computer program, the wing camber surface shape is found by 
integrating the surface slopes for all chordwise elements at 
selected span stations. Section values of drag {ACp times slope) , 
lift and pitching moment are then integrated spanwise to obtain 
wing force coefficients. 


Optimum comb in ation of loadings. - Equation 27, which defines the 
surface slope to support a specified lifting distribution, 
requires a description of the design pressure distribution. This 
is obtained from Lagrange* s method of undetermined multipliers, 
which provides a means of selecting an optimum combination of 
component loadings, yielding minimum drag subject to various 
imposed constraints on the allowable pressure distribution. 


A number of design point options are provided in the wing design 
program, involving the selection of loadings to be combined and 
the constraints desired. (The options are controlled by input 
codes, as described in the user's manual, part 2). 

The usual design conditions specified are: 

• Drag-due-ro- lif t of the wing be minimized at a 
given design lift coefficient, subject to an 
optional pitching moment constraint. 

• Effects of fuselage and nacelles be included in the 
design solution. 

• Constraints be applied to the design pressure 
distribution and/or local camber surface 
shape to provide physical realism. . 

The design point options are discussed on page 50. 

The optimum loading selection is an extension of the methods used 
in references 9 and 11, through the addition of the configuration- 
dependent loadings (due to fuselage and nacelles) and the addition 
of the pressure and ordinate constraint formulation. Also, more 
component loadings were added to provide flexibility in 
rearranging pressure distributions at the design point in order to 
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meet: the pressure constraint condition. The maximum number of 
leadings is 17, defined in the table on page 45 • 


The wing design method actually consists of four solution steps: 

• Flat wing solution 

• Calculation of aerodynamic characteristics of 

selected component loadings, 

• Optimization of loading combination 

• Camber shape calculation 


Flat wing so l utio n- - The first step in the wing design process is 
to obtain the flat wing solution for the given wing planform and 
Mach number, using the analysis form of the lifting pressure 
calculation described on page 5b . This solution is obtained to 
locate the wing aerodynamic center for use in computing the 
pitching moment at zero lift for the individual loadings; and 
also, because its drag-due-to-lift factor is an often-used 
reference point in judging twisted and cambered wing designs. 


Componen t l oadings characteri st ics. - For each selected component 
loading, section lift, drag, and pitching moment coefficients are 
evaluated at a series of span stations and integrated spanwise to 
obtain overall lift, drag, and pitching moment. In addition, the 
interference drag coefficients between loadings are computed, 
where the drag interference between loadings i and j is given by 


: D,ij C D,Ji 


■i- V AC .(L*,N*)(P-) ( L* ,N* ) A(L»,N») B(L»,N») 
pS AJ P,1 dX 


(28) 


+ T! AC )(!§•) (L*,N») A(L*,H«) B(L*,N») 

pS JL~t ox ^ 

where the summations are carried out over the wing planform. For 
the component loadings defined over the planform by analytical 
expressions, the loading is scaled by the appropriate power of 
either semispan or average wing chord, so that the resultant 
component loading lift coefficient is approximately one. 


Interference drag coefficients associated with the configuration- 
dependent loadings are special cases. The configuration- 
dependent loadings are those produced by fuselage upwash, fuselage 
volume asymmetry, and nacelles, as described on page 57 . They are 
calculated at discrete points over the wing planform, and 
interpolated linearly for a complete definition where needed . 
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TABLE / 

DESCRIPTION OF WING LOADING TERMS 

Loading 

Number Definition 


1. 

2 . 

3. 

4. 

5. 

6 . 

7. 

8 . 
9. 

10 . 

11 . 

12 . 

13 . 

14. 

15 . 

16. 
17. 


Uniform 

Proportional to x, the distance from the leading edge 

Proportional to y, the distance from the wing centerline 
2 

Proportional to y 

o 

Proportional to x 

Proportional to x(c - x), where c is local chord 
Proportional to x (1.5 c -x) 

Proportional to 2 (1 + 15 
Proportional to (1.05 c -x)®’® 

Elliptical spanwise, proportional to 
Proportional to x, the distance from the leading edge of an 
arbitrarily defined region 

A camber-induced loading proportional to the body bouyancy 
loading 

A camber-induced loading proportional to the body upwash loading 
A camber-induced loading proportional to the nacelle buoyancy 
loading 

The body bouyancy loading 
The body upwash loading 
The nacelle buoyancy loading 


'0 - y/|) 
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Since there is no camber surface associated with the 
configuration-dependent loadings, the interference drag 
coefficient between the "regular" loadings and the configuration 
dependent loadings has only one term. To ease the computational 
task of evaluating component and interference characteristics, the 
configuration-dependent loadings are assigned a corresponding wing 
slope definition of zero. 


The nacelle configuration-dependent loading has yet another 
difference, since the flow fields of the "regular” loadings act on 
the fixed geometry of the nacelles, producing an appreciable axial 
force. The lifting pressure acting on the wing lower surface is 
assumed to propagate along Mach lines from the wing to the nacelle 
and -*-he associated buoyancy force is computed. 


If 2 constraints are used, the camber surface 2 values of all 
loadings are interpolated at the planform locations where 2 
constraints are to be applied. These are used to constrain the 
optiirum loading combination by requiring the specified 2 value to 
occur at the constraint locations. 


If a fuselage is present, the wing solution proceeds outboard from 
the side-of-f uselage station, and carry-over lift and fuselage 
camber is used to define the force coefficients of the wing 
inboard of the side-of -f us elage. 

Fuselage contribut io n a nd "carry-over " li ft. - If a fuselage is 
present in the solution, the lifting pressure distribution of the 
component loadings inboard of the side-of-fuselage station is of 
the "carry-over" type, rather than the analytically defined 
loading. I.e., the loading on the wing area occupied by the 
fuselage is a dependent function of the loading defined outboard 
of the side-of-fuselage station. 

The carry-over pressure distributions are calculated for each 
component loading (excluding the configuration - dependent 
loadings) using the analysis form of the lifting pressure solution 
described on page 5^,. The carry-over pressure distributions are 
applied to the fuselage camberline to obtain the corresponding 
drag terms. 

In addition, the fuselage contribution to configuration lift, 
drag, and pitching moment is included in the design point force 
coefficients. The fuselage coefficients are calculated in the 
lift analysis program for the input fuselage geometry as described 
on page 62 and transferred to the wing design program. 

Optimization of loading combination. - Given the component- 

loadings and their interference drag terms, the aerodynamic 
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characteristics of a wing supporting any combination of these 
component loadings can be calculated in terms of load strength 
factors A | assigned to each. The total lift coefficient for n 
loadings is 


i=n 


C L “ 


C L,i A i 


i=l 


( 29 ) 


where 9^ j denotes the loading and A; its load strength factor. 
The total pitching moment coefficient at zero lift is C mo 


where 


C 
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i=n 



C 

mo ,i 
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and the total drag coefficient is 


i=n 


i=n 
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( 30 ) 


( 31 ) 


c 

where the terms DWN,i are the axial nacelle force coefficients 
for the component loadings. The contribution of each loading to 
the lifting pressure coefficient, the wing upper surface 
longitudinal pressure gradient and to the camber ordinate at point 
j can be summed to give 



i= 11 

E 


i=l 


C 




( 32 ) 


The load strength factors of the configuration-dependent loadings 
are 1.0, since those loadings are not variable. 


Hence, from the component loadings data, the drag, lift, and 
pitching moment characteristics of an infinite number of wing 
designs may be immediately calculated in terms of the A^ factors. 
The formal optimization of the wing loading involves the 
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specification of the set of A^’s which gives least drag subject to 
the imposed constraints. 


According to the Lagrangian method, the solution for the optimum 
A. values involves the system of linear equations shown in figure 
4.4-4, where matrix notation is used for compactness. The term 
* c P,i is the allowable lifting pressure coefficient corresponding 
to the i th pressure limiting constraint on the wing upper surface 
pressure. 


Partitions in the square matrix of the solution have been 
identified by number and letter in figure 4.4-4. The left, 
uppermost partition of this matrix (zone Al) and the adjacent row 
and column correspond to Grant* s original design solution 
(reference 14). The next m rows (zones C1-C7) and columns (zones 
A3-G3) are introduced by m local constraints on the lifting 
pressure coefficient, which are added (if necessary) in applying 
the pressure level constraint, as described on page 50 , The next 
q rows (zones D1-D7) and columns (zones A4-G4) are introduced by 
q constraints on the longitudinal gradient of wing upper surface 
pressure coefficient; these terms are added in applying the 
pressure gradient limiting condition. 


The next three rows (zones E1-E7) and columns (zones A5-G5) 
correspond to the constraints which set the loading factors to 
1.0 for the body bouyancy, body upwash, and nacelle buoyancy 
loadings. Each of these three rows and columns is present only if 
the corresponding configuration dependent loading is used. zones 
F1-F7 and A6-F6 correspond to r cpnstraints on ordinate, and 
appear only if ordinate constraints are actually specified by the 
user. The right-hand column (zones A7-G7) and bottom row (zones 
G1-G6) correspond to the constraint on pitching moment coefficient 
at zero lift; this constraint is also optional. 


The unknowns of the design solution are the n values of the 
loading factors and the (m+q+r+5) Lagrange multipliers ( Ai ) • 
These are contained in the left-hand-side column matrix. 


The seven partitions of the right-hand- side column matrix contain; 
either the negative sum of the nacelle interference drag produced 
by wing lift and fuselage carry-over drag, or zero if not used; 
the design lift coefficient less fuselage lift coefficient; the 
pressure level constraint values applied to the lifting pressure 
coefficient (these are generated automatically by the design 
computer program, if they are necessary) ; the pressure gradient 
constraint values (also applied automatically if they are 
necessary) ; the values of for the three configuration-dependent 
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loadings; the constrained values of camber-line ordinate; and the 
design pitching moment coefficient at zero lift less fuselage 
pitching moment at zero lift. 


The maximum number of equations in figure 4.4-4 is 34. This 
number of equations can be reached if 17 loadings are combined, a 
pitching moment constraint is used, and the maximum number of 
constraints on lifting pressure coefficient are required. 


Design camber surface. - Given the set of A, factors, the design 

pressure distribution is known and the resulting camber surface is 
calculated using equation 27 on page 40 


Design_goint_gptiqns. - Due to options in the wing design program 
tc limit the solution extent (e. g. , no C mo constraint or no 
ordinate constraints or no pressure constraints) , some parts of 
the solution matrix on page 49 may not be used. However, for 
purposes of describing the solution options, it is assumed that 
all options are chosen. In that case, the wing loading solution 
proceeds as follows, using repeated applications of the 
corresponding matrix: 

(1) The wing loading having minimum lift-dependent drag with 
a constraint on design lift coefficient is defined. 
This solution will include ordinate constraints if any 
have been specified. 

(2) A family of wing loadings having minimum lift- dependent 
drag with a constraint on design lift coefficient and a 
series of constraints on pitching moment coefficient at 
zero lift are defined. This series of solutions is 
presented in the form of drag due to lift factor. Kg, 
versus C mo , and is referred to as the "bucket” plot. 
Ordinate constraints are not included. 

(3) The wing loading having minimum lift-dependent drag with 
constraints on design lift coefficient wing upper 
surface pressure and camberline ordinates is found. 
This solution may require a cyclic process for the 
pressure constraints that begins with (1) above. If (1) 
satisfies the pressure criteria on both gradient and 
level everywhere on the planform, then this loading is 
set equal to (1) . If (1) does not satisfy the pressure 
limits everywhere, then a pressure constraint is applied 
at the wing planform location where the pressure 
gradient is most strongly violated, and a new solution 
for wing loading is found. Its pressure distribution is 
tested against the gradient limit, and another 
constraint is imposed if the limit is violated. Each 
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cycle adds another pressure gradient constraint to the 
solution in addition to the lift coefficient constraint. 
When pressure gradient is satisfied, then pressure level 
is similarily tested against its limit. Constraints on 
pressure level are also added cyclically until solution 
pressure level satisfies the level limit everywhere. 
Cycling continues until either the solution load 

distribution everywhere satisfies the pressure gradient 
and level limit or until the number of solution 
constraints equals the maximum allowable number. In the 
latter case, the limit level for pressure gradient is 
increased by 20%, and the cycle is started again. 

It has been found that pressure gradient and pressure 
level constraints imposed early in this cyclic solution 
can subsequently become unnecessary as other pressure 
constraints are imposed. Constraints unnecessary for 
satisfying pressure gradient and level criteria penalize 
the wing design unnecessarily by increasing drag. Such 
constraints can be detected by the algebraic signs of 
their solution Lagrange multipliers. Accordingly, a 

test for unnecessary pressure constraints has been added 
to the solution cycle. Such constraints are removed 
when detected and cycling continues. 

<4) The wing loading having minimum lift -dependent drag with 
constraints on design lift coefficient, camberline 
ordinates, and zero-lift pitching moment coefficient is 
defined. 

(5) The wing loading having minimum lift-dependent drag 

with constraints on design lift coefficient, pitching 
moment coefficient at zero lift, camberline ordinates, 
and wing upper surface pressure is found. The latter 
type of constraints are imposed, if necessary, in the 
cyclic fashion of (3). 

Examples of loading solutions (1) , (2) , (4) , and (5) are shown in 
figure 4.4-5. For this case, loading solution (1) satisfies the 
wing upper surface pressure constraint, and is therefore identical 
to leading solution (3) . 

The effect of the number of loadings on the "bucket 11 plot and on 
loading solution (5) is illustrated in figure 4.4-6. Increasing 
the number of loadings from 3 to 11 results in a substantial 
theoretical drag decrease; however, this decrease is reduced 
substantially when pressure gradient constraints are added. 
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Wing upper surface C constraint: C >0.7 C 

■ J ^ p vacuum 



<2 

mo 
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FIGURE 4.4-5. ~ TYPES OF L OA DING SOL U TIONS 




NOTE: 

Flagged symbols denote 
solutions satisfying 
dCp 

< 0.0259* 
dX 

on the wing upper surface, 
including thickness Cp 


0.6 


Flat wing 


0.5 


Families of solutions 
without pressure 
level or gradient 
constraints: 



*Maximum for 3 loading solution, 
including thickness pressures. 


FIGURE 4.4-6-WING DESIGN STUDY: CONSTRAINTS ON BOTH C mo AND UPPER SURFACE 
PRESSURE GRADIENT 
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A discussion of the loading selection, constraint application, and 
the resultant camber surface optimization is presented in part 2 
(user's manual). 


Lift Analysis Case 


Loading for a gi ven camber surface. - Since the R function for the 
field point element is zero (the element generates no downwash or 
upwash upon itself) , equation 27 can be rearranged to solve for 
the lifting pressure coefficient in terms of the field point slope 
and upstream influences: 

» dz 

&c p (L»,N*) 

+ Y ~ L > N * ~ N ) A(L,N) B(L,N) C(L,H) AC : (L,K) 

T 

Therefore, the lifting pressure distribution, ACp, 
determined for a wing of arbitrary surface shape provided the 
calculations are performed in the proper sequence. The order of 
calculation of ACp(L*,N*) is from apex aft, and from the 
centerline to the right hand wing tip for each L* row. In this 
fashion, all pressure coefficients within the Mach forecone from 
any element will have been previously obtained and no unknown 
pressure coefficients arise in the summation. The element 
fractions A, B, and C are as defined for equation 27. 


(33) 

can be 


Theoretically, ACp(L*, N*) defined by equation 33 is the pressure 
coefficient at the aft mid- point of the L* , N* element. The 
average pressure coefficient for the element, needed in subsequent 
calculations, is calculated by one of two alternative methods: 

1) An approximate average pressure coefficient is 
interpolated from the aft mid- point Acp value and the 
ACp value of the element immediately ahead. (No 

interpolation is performed if L*, N* is a leading edge 
element) . This is the method used in the original 
formulation of the computer program. However, it was 
found that osillations in the ACp values occurred near 
the wing leading edge. So, a 9 point smoothing equation 
was applied after unsmoothed pressures were computed for 
the entire wing, which essentially removed the 

oscillations and brought the smoothed pressure 

distribution into good agreement with established 
results from other theoretical methods. The smoothing 
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equation is of the form 

( 3*0 

0.2A{L # -4)AC ) + 0.4a(L»-3)AC (L*-3) + 0.6A(L*-2)AC (L*-2) 

P P P 

♦ 0.8A(L«-l)AC p (L«-l) + A(L»)AC p (L*) + 0.8A(L» +l)AC p (L* + l) 

_ + 0 . 6A ( L # +2 ) AC p ( L*+2 ) + 0.1*A(L # +3)AC p (L*+3) + 0.2A(L*+4)AC p (L*+4) 

AC p (L*) “ 0.2A(L»-U) + 0.4 a{L # -3) + 0.6a(L*- 2) + 0.8A(L»-l) + A(L*) + 0.8A(L # +l) 

+ 0.6a(L»+2) + O.UA(L # +3) + 0.2A(L*+4) 


where the A. values account for element fractions. A is 
zero if the corresponding element is ahead of the wing 
leading edge, or aft of trailing edge in the case of a 
subsonic trailing edge. The supersonic trailing edge 
solution is treated as a special case, and the trailing 
edge is extended four elements to provide pressure 
coefficient data to fill out the smoothing equation. 


2) The other average ACp method is an improved technique, 
which uses an aft element sensing approach. This method 
involves solving for preliminary ACp results for a given 
L* , N* element and the element immediately aft, then 
following it up with a second calculation to refine the 
preliminary results. The procedure is detailed in 
reference 12, but is summarized below and in figure 4.4- 
7. 


a) 


c) 


Calculate preliminary Acp values for a given L* = 
constant row. Designate as ACp^ a (L*,N*) 


Calculate preliminary Acp values for the element row 
immediately aft - L*+1, using Ac values. 

Designate as AC (L*,N*). P ' a 

P » D 


Finalize Acp values for the original L* = constant 
row from one of the two following equations. 

For leading edge elements: 


AC (L*,H») • 



AC (L»,N») 
p,a 


. JL f A(L«.K») 

2 [l + A(L»,N«) 


AC . (L»,N») 

P,D 


(35) 


For all other elements: 




AC (L*,N*) 
P»a 


AC ( L* , N* ) 

p,b 


( 36 ) 


The aft element sensing technique produces a substantial reduction 
in pressure coefficient oscillation- In cambered wing test cases 
run with the method, however, some oscillation was found to 
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final AC for 
all elements 
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-=l_l 
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& 
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I 

“? 1 — Preliminary Ac„ for (L*, N*), AC„ „ (L*, N*> 

■ * ■ r P ; a 




N — Fairing of preliminary results 
' — Integration of preliminary AC^ 


J_Z i i i 



FIGURE 4.4-7.— APPLICATION OF AFT-ELEMENT SMOOTHING TECHNIQUE (FROM REF. 12) 



persist. Therefore, a 3 term smoothing equation was added (after 
all wing pressures are calculated) , of the form: 

0.5A(L»-l)AC ( L*-l ) + AC (L* ) + 0. 5A(L*+1 )AC (L*+l) 

AC p ( Lll )= 0.5A(L*-1) + 1.0 + 0 . 5A ( L # +l ) “ 

In the computer program one of the two alternate smoothing 
methods is selected by an input code. Both produce essentially 
the same answers • When the pressure limiting option is used 
(discussed later) , the aft-sensing smoothing technique is 
automatically selected. 


For the analysis solution, the pressure coefficients for all 
elements must be calculated. The force coefficients are, 
therefore, calculated from direct summations of local pressures 
applied to each element, rather than employing a spanwise 
integration as in the design case. Lengthwise and spanwise lift 
distributions are obtained by summing the lift in the 
corresponding element rows. 


No " leading-edge-suction" force is included in the drag summation, 
which accounts only for the lifting pressure coefficients acting 
on the local wing slopes. 

Analysis of complete configuration. - The preceding discussion 

describes the theory used in the lifting pressure calculations for 
a wing camber surface at a selected angle of attack. The analysis 
pregram may also be used to calculate lifting pressure 
distributions and force coefficients for complete configurations 
over a range of angles of attack, adding in the effects of 
fuselage, nacelle, canard and/or horizontal tail, as applicable. 


The program actually carries two solutions along: one for the 
configuration at its input angle of attack, the other the 
incremental solution per degree angle of attack (called the flat 
wing solution) . The interference terms associated with the two 
solutions acting on the other surface shape are also calculated. 
The summation of these effects into the drag polar and other force 
coefficients is performed by superposition, as described on page 
64 


Calculation of the complete configuration lifting pressure 
solution involves up to 7 principal tasks: 


Isolated fuselage upwash field 

Nacelle pressure field acting on wing 

Pressure field due to asymmetrical fuselage volume 
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• Wing/canard solution in presence of fuselage upwash 
field- Effects of the wing pressure field acting 
on nacelles, are also calculated 

• Fuselage lift distribution in presence of wing 
downwash field 

• Horizontal tail solution in presence of fuselage 

and wing flow fields, 

• Superposition of solutions 

The calculation logic of these tasks is described in the following 
sections : 


Isolated f us elage upwash field. - The fuselage upwash (or 

downwash) field in the plane of the wing, canard or horizontal 
tail is obtained from slender body theory. From reference 15, the 
velocity potential of a slender body of revolution, shown in 
figure 4.4-8, is: 

/ x 

dA(x) / s . r . / Mr, 2 / \ cos 9 

-as— £n(x - x i )dx i + Ca+e(x)lR (x) — — 


The ncndimensional radial velocity component is 


i_ dA{x) _ a + e(x)_ R 2 ( x ) cos g 


U r 


2 ttt dx 


(39) 


and the tangential component is 


1 * - -[a ♦ c(x)]-^M sin 0 


(HO) 


The vertical velocity component (upwash angle) is 


H»r cos 9 -4 +e sin 6 


(41) 


Substituting for and 
gives 

w = _1 dA(x) z_ _ [ + , Rllx) 

U 2n dx 2 1 WJ 2 


and converting to cartesian coordinates 

2 , 


2 2 
z - y 


(h2) 


When the fuselage is at input incidence, the fuselage upwash is 

(43) 


_ 1 dA(x) z < , R 2 (x) / z 2 - y 2 
U o 2 tt dx r 2 " ' ; “^2 \ r 2 


and the incremental upwash per degree angle of attack increase is 


* _ d w/u __ R^(x) / z 2 - 

U ~ dot 2 1 ~2 


(kb) 
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Mach 'S 
line 


Upwash propagates along Mach lines. 
Upwash at arbitrary point P is: 



w 

u 


± 

27 r 


dA(x) z_ , . R 2 (x) 

dx r 2 r 2 



= _ R 2 (x) / z 2 - y 2 \ 
r 2 V r 2 / 


Where: 

u Q = upwash at input fuselage incidence (0:= 0) 

u = upwash increment per unit a 

A(x), R(x) = Fuselage cross-sectional area, radius 

e (x) = Fuselage incidence due to camber 
{relative to reference plane) 

z,y,r = location of point at which 
upwash is calculated 

a = angle of attack of reference plane 
U = free stream velocity 


FIGURE 4.4-8.-FUSELAGE UPWASH CALCULA T/ON 
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Slender body theory assumes that the upwash field propagates 
normal to the body centerline. In the analysis program, the 
fuselage upwash is propagated along Mach lines, as required by the 
characteristics of supersonic flow. 


Nacel le pressure fiel d acting o n wing. - The pressure field caused 
by the nacelles on the wing is calculated by the Whitham solution, 
as described in Section 4.3. The nacelles are assumed to be 
bodies of revolution, and the pressure signature due to each 
nacelle is calculated for a series of spanwise stations from wing 
root to tip. The composite pressure signature at a given spanwise 
station is the sum of the individual nacelle signatures, with all 
pressure coefficients doubled to account for reflection from the 
wing. Either "wrap” or "glance" nacelles may be used, as 
discussed in Section 4.3. 


Press ure fiel d due to asymmetric fu s ela ge volu me. - Another term 

in the superposition approach is associated with the growth in 
fuselage area above and below the wing. If the area growth is 
asymmetric (e.g., a low wing configuration), a differential 
pressure across the wing plane is created which is additive to the 
effect created by the fuselage upwash field. 


An approximate method of computing the asymmetric fuselage 
pressure field is used, using the area distributions shown in 
figure 4.4-9. To the actual forebody area distribution is added 
the growth in fuselage area in the wing region relative to the 
area at the leading edge of the wing-f uselage intersection. This 
is dene for both the below- amd above-wing area distributions. 
Pressure signatures for both body representations are computed 
using the Whitham technique at a series of spanwise stations, 
doubled to account for reflection from the wing, and differenced 
to get the lifting pressure distribution due to asymmetric 
fuselage volume. 


If the wing trailing edge is subsonic, the pressure fields above 
or below the wing revert to that corresponding to the actual 
fuselage area distribution in the region aft of a Mach line from 
the trailing edge of the wing-fuselage intersection. 


Use of the asymmetric fuselage calculation is controlled by an 
input code in the computer program. If the asymmetric calculation 
is net requested, the fuselage pressure field corresponding to a 
mid-wing arrangement is calculated, so that there will always be 
a thickness pressure field due to the fuselage, if present, for 
use in limiting pressure calculations (discussed later). 
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Wing/Canard Solution in Presence of Fuselage and Nacelles 


The wing or canard lifting pressure solution in the presence of 
fuselage and/or nacelles is performed as described for the 

isolated wing, with the following alterations: 

• The local surface angles of attack are increased by the 

fuselage upwash values for the purposes of computing 
pressure coefficients. (To compute drag, the pressure 
coefficient is applied to the local surface slopes, 
only) . 

• In the region of the wing covered by the fuselage, the 

wing slopes are zeroed. This reflects the fact that the 
wing lift distribution in this region is of the "carry- 
over" type, only. (To compute drag, the carry-over lift 
is applied to the fuselage camber line shape) 

• The canard pressure distributions are computed in the 

same fashion as the wing pressure distribution, with the 
canard also creating a downwash (or upwash) field for 
the wing. 


The nacelle or asymmetric fuselage volume pressure fields are 
superimposed upon the wing elements to obtain the effect of the 
nacelles or fuselage asymmetry on the wing- If nacelles are 
present, the buoyancy drag of the wing lifting pressures acting on 
the nacelle area distribution is also calculated. This term is 
obtained by transferring local wing pressures aft along Mach lines 
to act upon the nacelle area growth, in a fashion paralleling the 
wing-thickness-on- nacelle term of the near- field wave drag 
program. 


Fuselage l if t distribution in p re sence of wi ng d ownwash fi eld . - 
From slender body theory , the fuselage lift distribution is given 
by 


eg • R(x) = — ^ [a + e(x) ] + A(x) _JL.[e(x) ] 

dx 


(U5) 


using the notation of figure 4.4-8. 

The first term of the equation. is the slender body term associated 
with a straight body at angle of attack; the second term is due to 
the curvature change associated with body camber. 
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In computing the fuselage lift distribution, the downwash (or 
upwash) from the wing/canard is included in the local angle of 
attack term, e(x). For drag computations, the lift is applied only 
to the physical angle of attack value, a+e(x). 


Fuselage force coefficients (lift, drag, and pitching moment) are 
calculated from the slender body lift distribution and converted 
to a wing reference area basis for summation into the complete 
configuration characteristics. 


Horizontal tail solu ti on. - The horizontal tail lifting pressure 
distribution is calculated in the presence of the fuselage upwash 
field and wing/canard downwash field. The solution employs the 
same logic of partitioning the tail into exposed and carry-over 
regions used in the wing and canard pressure computations. 


Since the horizontal tail is assumed to be aft of the wing, the 
fuselage and wing/canard upwash (or downwash) field is first 
calculated, then the tail lift distribution and force coefficients 
are computed for the desired tail angles of attack. These are 
summed into the wing/fuselage, etc., solution for each tail angle, 
resulting in a set of force coefficients for each tail setting in 
the same fashion as wind tunnel data are obtained. 


Calculation of the wing/canard downwash field is performed by 
extending the wing grid system aft to include the horizontal tail. 
With all pressure coefficients on the wing or canard previously 
established, the effect of the wing/canard on the horizontal tail 
is obtained from 

Wing 

A u p = Z) R (L* - L, N* - N) A (L, II) B(L, N) C(L, N) AC (L, N ) 

T P 

Canard 

+ £ R (L* - L, N*- N ) A(L, N) B(L, N) C(L, N ) AC (L, N) 

T P 

in which the wing, canard, and tail are all assumed to be located 
in essentially the same plane. 


Calculation of the horizontal tail lifting pressures are then 
performed using the analysis form of the lifting pressure equation 


ACp «■*. ■*> = 'I E »*) 

(bj) 

Tail _ 

+ i y* R (L* - L, N* - N) A(L, II) B(L, N) C(L, N) AC (L, N) + A, T _ 

IT ““ P W t L> 
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with the wing/canard term added to the summation of upstream 
effects. As in the wing/canard pressure solution, the fuselage 
upwash is added to the physical tail slope for the purposes of 
computing lifting pressures, but removed for the calculation of 
the drag term. The carry-over lift of the horizontal tail is 
applied to the fuselage camber line to obtain the drag 
interference of the tail on the fuselage. 


The effect of wing downwash on the horizontal tail lifting 
pressures is significant. The data of figure 4.4-10 show the 
theoretical drag and pitching moment increments at constant total 
lift for a typical configuration, with and without wing downwash 
included. (Comparisons of the theoretical calculations with 
experimental data are presented in the user’s manual, part 2.) 

C anar d or horiz o nta l tail downwash shift. - The lift analysis 

program contains a feature to shift the downwash field of wing or 
canard in the horizontal plane. The basic theoretical solution 
assumes the downwash field moves directly aft, instead of 
following the fuselage contour, as illustrated in figure 4.4-11. 
The shift feature can have significant influence on the calculated 
interference characteristics, as shown in figure 4.4-12. Program 
options allow for: 

• No shift 

• Shift according to the side -of- fuselage span stations of 
wing, canard, or tail (normal option). 

» Shift according to an input dimension 


Superposition of Solutions 


The arrangement of terms in the superposition of the complete 
configuration force coefficients is: 

C T = 2 C T + a £ AC (48) 

L L a = 0 L 


c = £c 


m 


ra 


+ a £ AC 


a = 0 


m 


(U9) 


c. = 2c 


+ K„ 


a - 0 


(C L - 


) + K, 


a = 0 


(C L- 


= 0 


(50) 
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M = 2.7, C L = 0.10 


Coefficients based on wing area 
Increments from tail on/tail off data 






Flow angle at tail 
downwash plane, deg 



FLOW ANGLE AT 
TYPICAL DOWNWASH PLANE 


Wing side of 
fuselage 


Upwash station 



station 


FIGURE 4.4-11 —DOWNWASH PATH CORRECTION 



LEGEND 



FIGURE 4.4-12 -EFFECT OF CANARD DOW NW ASH SHIFT ~ M = 2.7 
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where the lift, drag, and pitching moment coefficients at«= 0 
correspond to the individual fuselage, wing, tail, etc., solutions 
at input configuration incidence, and the AC^ and AC m 
coefficients correspond to the incremental (flat wing) solutions 
per degree angle of attack. 


The drag coefficient equation includes the interference terms 
between the solution at input incidence and the incremental 
solution per degree. With subscripts C and F denoting input 
incidence and the flat wing terms, respectively, the drag 
coefficient equation is: 


or 


C D = % + (C D poc 


+ CL. ) — + — C_ a 
D COF a F a F L F 


(C 


+ C 


C = C 
D D, 


where 


FOC 


) 


COF 


,c - - v * 


0 . 017 U 5 (C T - C T )' 

L L c 


**F "F 

oc = configuration angle of attack 


(51) 

( 52 ) 


a = flat wing incidence, .01745 radian 

r 


-l f 


■ s 2 Ac p c 
= A 


D c = i £AC P 0 A E “C 


D, 


= “F \ 


with a e = local element area corresponding to ac^ 
S = reference area 


0 ^ = local element slope 


2 interference drag of flat wing 

C D = lf£ AC p A E a C pressures acting on cambered 
F0C F wing slopes, per degree. 


: °C0F = “ F % 


interference drag of cambered 
wing pressures acting on flat 
wing incidence, per degree. 


with the summations carried out over the wing planform. 


Therefore, for each of the solutions involved in the 

superposition, it is necessary to calculate both an input 
incidence solution, a flat wing solution, and the interference 
drag terms between the two solutions. All of the force 

coefficients are referred to an input area, moment center, and 
moment reference length in the summations. 
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Limiti ng p res sure solution . - Linear theory imposes no 

restrictions on the allowable wing pressure coefficients. 
Particularly in the case of subsonic leading edge wings, large 
upper surface pressure coefficients may be computed near the 
outboard wing leading edge that exceed (are more negative than) 
pressure coefficients corresponding to a vacuum. These pressure 
coefficients can occur at moderate wing angles of attack, due to 
the strength of the upwash generated by the inboard wing. 


Several investigators have examined experimental data to define 
the minimum attainable upper surface pressure coefficients. One 
such correlation is shown in figure 4.4-13 (from reference 16), 
which indicates a physical pressure coefficient limit of 
approximately .8 vacuum. The mechanism of this limiting is 
associated with local leading edge flow separation. However, the 
phenomenon may be approximated in the wing analysis program by 
limiting the calculated linear theory pressure coefficients to a 
specified fraction of vacuum Cp. 


In the analysis program, pressure limiting is an optional feature, 
controlled by an input code. If limiting is requested, then a set 
of configuration angles of attack is required for the solution, 
since superposition will not apply after limiting occurs. Also, 
a definition of the wing thickness pressures is required 
(transferred over from the near field wave drag program module) , 
since it is the sum of lifting plus thickness pressures that is 
limited, rather than the lifting pressure alone. 


Pressure coefficients, as calculated, are separated into upper and 
lower surface values. Thickness pressures (wing plus fuselage) 
are added and the upper surface pressure coefficient tested 
against the limit pressure coefficient. If the limit is exceeded, 
the calculated lifting pressure coefficient is adjusted to a level 
which causes the total upper surface pressure coefficient to natch 
the limit value. 


Comparisons of the "limiter" calculations with measured force and 
pressure coefficient data are presented in figures 4.4-14 through 
4.4-18 for 2 arrow wings of references 1 and 2. Both were 70° 
sweptback wings, having 3 percent biconvex airfoils, and wing 
design lift coefficients of 0 (flat) and .08. 


The basic thickness pressure comparison is shown in figure 4.4-14. 
Comparisons of upper and lower surface pressure coefficients at 
lifting conditions of Cl% .13 and .25 are shown in figures 4.4-15 
and 4.4-16. In the case of the flat wing (figure 4.4-15) , several 


69 



fractions of vacuum pressure coefficient are illustrated; the 
design C-^ = .08 wing data are shown for .7 vacuum limiting only. 


The limiter feature exhibits considerable improvement, in the 
detail pressure coefficient comparisons over unlimited linear 
theory as Cj- increases. 

Force coefficient comparisons for the two wings are shown in 
figures 4.4-17 and 4.4-18. Pressure limiting at approximately .7 
vacuum improves the comparisons appreciably at the higher 
values. 
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FIGURE 44-15. —PRESSURE COEFFICIENT COMPARISON, WING 1 (FLAT WING) 
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FIGURE 4.4- 16. -PRESSURE COEFFICIENT COMPARISON, WING 2 

(TWISTED AND CAMBERED WING, DESIGN C, = .08) 
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5.0 INTERACTIVE GRAPHICS 


Interactive graphics for use with the design and analysis system 
are optional, employing the NASA-LRC cathode ray tube (CRT) 
display and associated software. Use of the graphics option is 
requested by an executive control card (described in part 2, 
user's manual) . 

The principal uses of the graphics routines in the design and 
analysis system are to display the configuration, edit input 
geometry, and to display and/or alter the basic program 
calculations. Limited capability to redirect the system 
calculation sequence is available from the CRT console during 
program execution. 

The CRT display and program coding for the interactive graphics 
setup are based on the NASA-LRC system. However, all display 
portions of the coding are subroutined or overlaid from the basic 
programs, so that the system could be readily converted to other 
CRT arrangements. 

Details of the interactive graphics portion of the design and 
analysis system are given in part 2 (user's manual). 
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